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Outline  of  the  talk 


1 .  Some  (well  known)  Lie  group  methods  for  linear  problems  (Fer  and 
Magnus  expansions). 


2.  Schemes  based  on  triangular  matrices  (splitting  +  solvability). 


3.  Some  methods  and  practical  issues  in  their  construction 
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1  Lie  group  methods  (linear  probiems) 


Let  us  consider  a  linear  matrix  ODE  evolving  in  a  Lie  group  Q 


with  A  :  [ro,°°[x  ^  g  smooth  enough, 

g:  Lie  algebra  associated  with  ^ 

Examples  of  SL(n),  0(n),  SU(n),  Sp(n),  SO{n) 


7(r)  G  Lie  group  ^  if  A(r)  g  Lie  algebra  g 

*  There  are  several  schemes  preserving  this  feature  (Magnus,  Per, 
Cayley,...) 
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1 .1  Magnus  expansion 


For  the  equation 


Y'=A(t)Y,  Y(to)  =  I. 


*  Magnus  (1 954)  proposed 


oo 


k=l 


(1) 
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1 .1  Magnus  expansion 


For  the  equation 


Y'=A(t)Y,  Y(to)  =  I. 


*  Magnus  (1 954)  proposed 


oo 


k=l 


(1) 


with  log(y(r))  satisfying 


oo 


Q.'  =  d exp A{t)  =  £  l^adJiA(r) 

k=0 


=  0: 


(2) 
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1.1  Magnus  expansion  (II) 


ad^A  =  A 

ad^A  =  [n,ad^^A] 


[^^,A]  =  QA  -AQ. 


and  Bk  are  Bernoulli  numbers. 
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1.1  Magnus  expansion  (III) 


First  terms  in  the  expansion  (A/  =  A(r/)): 


^i(0  =  [  A{t\)dti 

Jto 

^2(0  =  o  /  /  <^^2[2il,A2] 

^  Jto  Jto 

^3(0  =  A  /  [2i2,A3]]  +  [A3,  [A2,Ai]]) 

o  Jto  Jto  Jto 


e^(0  g  g  even  if  the  series  Q.  is  truncated 
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1.1  Magnus  expansion  (III) 


First  terms  in  the  expansion  (A/  =  A(r/)): 


e^(0  g  g  even  if  the  series  £l  is  truncated 


*  Expansion  widely  used  in  Quantum  Mechanics,  NMR  spectroscopy, 
infrared  divergences  in  QED,  control  theory,... 
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1 .1  Magnus  expansion  (IV) 


Magnus  as  a  numerical  integration  method  (Iserles  &  Norsett,  1997) 
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1 .1  Magnus  expansion  (IV) 


Magnus  as  a  numerical  integration  method  (Iserles  &  Norsett,  1997) 


Two  critical  factors  in  the  computational  cost  of  the  resulting  algorithms: 
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1 .1  Magnus  expansion  (IV) 


Magnus  as  a  numerical  integration  method  (Iserles  &  Norsett,  1997) 


Two  critical  factors  in  the  computational  cost  of  the  resulting  algorithms: 
(1)  Evaluation  of  exp(^2) 

(Moler  &  Van  Loan,  Celledoni  &  Iserles,...) 
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1 .1  Magnus  expansion  (IV) 


Magnus  as  a  numerical  integration  method  (Iserles  &  Norsett,  1997) 


Two  critical  factors  in  the  computational  cost  of  the  resulting  algorithms: 

(1)  Evaluation  of  exp(^2) 

(Moler  &  Van  Loan,  Celledoni  &  Iserles,...) 

(2)  Number  of  commutators  involved  in  the  expansion 
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1 .1  Magnus  expansion  (IV) 


Magnus  as  a  numerical  integration  method  (Iserles  &  Norsett,  1997) 


Two  critical  factors  in  the  computational  cost  of  the  resulting  algorithms: 

(1)  Evaluation  of  exp(^2) 

(Moler  &  Van  Loan,  Celledoni  &  Iserles,...) 

(2)  Number  of  commutators  involved  in  the  expansion 

To  reduce  this  number  is  particularly  useful  the  concept  of  graded  free 
Lie  algebra  (Munthe-Kaas,  Owren  1999) 
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1 .1  Magnus  expansion  (V) 
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1 .1  Magnus  expansion  (V) 


*  Numerical  schemes  based  on  Magnus  up  to  order  8  have  been 
constructed  involving  the  minimum  number  of  commutators  in  terms  of 
quadratures  and/or  univariate  integrals. 
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1 .1  Magnus  expansion  (V) 


*  Numerical  schemes  based  on  Magnus  up  to  order  8  have  been 
constructed  involving  the  minimum  number  of  commutators  in  terms  of 
quadratures  and/or  univariate  integrals. 

*  Efficient  in  applications 
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1 .2  Other  schemes 


Fer  expansion.  Built  by  Francis  Fer  (1 958). 
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1 .2  Other  schemes 


Fer  expansion.  Built  by  Francis  Fer  (1 958). 


*  Referred  erroneously  in  the  (mathematical  physics)  literature  (e.g., 
Wilcox  1967),  but ... 
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1 .2  Other  schemes 


Fer  expansion.  Built  by  Francis  Fer  (1 958). 


*  Referred  erroneously  in  the  (mathematical  physics)  literature  (e.g., 
Wilcox  1967),  but ... 


*  proposed  (as  an  exercise!)  by  R.  Bellman,  ‘Introduction  to  Matrix 
Analysis’,  1960,  page  204: 


NEW  NUMERICAL  INTEGRATORS  BASED  ON  SOLVABILITY  AND  SPLITTING  -  20 


C ◄ ►DOX 


1 .2  Other  schemes 


Fer  expansion.  Built  by  Francis  Fer  (1958). 

*  Referred  erroneously  in  the  (mathematical  physics)  literature  (e.g., 
Wilcox  1967),  but ... 

*  proposed  (as  an  exercise!)  by  R.  Bellman,  ‘Introduction  to  Matrix 
Analysis’,  1960,  page  204: 

“The  solution  of  dX/dt  =  Q{t)X,  X(0)  =  /,  can  be  put  in  the  form 

•  •  • ,  where  P  =  /q  Q{s)ds,  and  Pn  =  /q  Qnds,  with 


The  infinite  product  converges  it  t  is  sufficiently  small.” 

(See  also  Mathematical  Reviews  21  2771,  review  done  by  R.  Bellman) 
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1 .2  Other  schemes  (II) 


*  used  as  an  (analytic)  procedure  in  perturbation  theory  in  Quantum 


Mechanics  by  Klarsfeld  &  Oteo  (1989),  but... 
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1 .2  Other  schemes  (II) 


*  used  as  an  (analytic)  procedure  in  perturbation  theory  in  Quantum 


Mechanics  by  Klarsfeld  &  Oteo  (1989),  but... 


*  numerical  integration  method  built  by  Iserles  (1984). 
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1 .2  Other  schemes  (II) 


*  used  as  an  (analytic)  procedure  in  perturbation  theory  in  Quantum 
Mechanics  by  Klarsfeld  &  Oteo  (1989),  but... 

*  numerical  integration  method  built  by  Iserles  (1984). 

*  This  class  of  methods  can  actually  be  built  from  Magnus. 
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1 .2  Other  schemes  (II) 


*  used  as  an  (analytic)  procedure  in  perturbation  theory  in  Quantum 
Mechanics  by  Klarsfeld  &  Oteo  (1989),  but... 

*  numerical  integration  method  built  by  Iserles  (1984). 

*  This  class  of  methods  can  actually  be  built  from  Magnus. 

*  They  require  the  computation  of  several  matrix  exponentials. 
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1 .3  Methods  based  on  the  Cayley  transform 


Let  us  suppose  that  Y'  =A{t)Y  is  defined  in  a  /-orthogonai  Lie  group, 


Oj{n)  =  {A  G  GL„(R)  :  A^JA  =  J} 


J\  constant  matrix 
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1 .3  Methods  based  on  the  Cayley  transform 


Let  us  suppose  that  Y'  =A{t)Y  is  defined  in  a  /-orthogonai  Lie  group, 


Oj{n)  =  {A  G  GL„(R)  :  A^JA  =  J} 


J\  constant  matrix 

Examples:  orthogonal  group  {J  =  I),  symplectic  group,  Lorentz  group 

(/  =  diag(l,-l,-l,-l)). 

Solution: 
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1 .3  Methods  based  on  the  Cayley  transform  (II) 


with  C{t)  G  oj{n)  satisfying  (Iserles  2001) 


dC 

dt 


t  >  to, 


C(to)  =  0. 


efficient  methods  without  matrix  exponentiais! 
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1 .3  Methods  based  on  the  Cayley  transform  (II) 


with  C{t)  G  oj{n)  satisfying  (Iserles  2001) 


dC 

dt 


t  >  to, 


C(fo)  =  0. 


efficient  methods  without  matrix  exponentials! 

In  fact,  one  can  also  combine  Magnus  with  Fade  to  avoid  the  use  of 
matrix  exponentials  in  /-orthogonal  groups! 
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1 .3  Methods  based  on  the  Cayley  transform  (II) 


with  C{t)  G  oj{n)  satisfying  (Iserles  2001) 


dC 

dt 


t  >  to, 


C(fo)  =  0. 


efficient  methods  without  matrix  exponentials! 

In  fact,  one  can  also  combine  Magnus  with  Fade  to  avoid  the  use  of 
matrix  exponentials  in  /-orthogonal  groups! 

*  It  is  possible  to  construct  methods  which  are  more  efficient  than  those 
based  on  the  Cayley  transform  (Blanes,  C.,  Ros  2002). 
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1.4  Summary 


*  These  methods  require  the  evaluation  of  one  or  severai  matrix 


exponentials 
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1.4  Summary 


*  These  methods  require  the  evaluation  of  one  or  several  matrix 
exponentials 


They  are  expensive  when  n  is  very  large 
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1.4  Summary 


*  These  methods  require  the  evaluation  of  one  or  several  matrix 
exponentials 


^  They  are  expensive  when  n  Is  very  large 

*  In  some  cases,  If  the  exponential  Is  approximated  by  rational  functions 
the  method  does  not  preserve  the  Lie  group  structure, 

in  particular,  when  Q  =  SL(n) 
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1.4  Summary 


*  These  methods  require  the  evaluation  of  one  or  several  matrix 
exponentials 


^  They  are  expensive  when  n  Is  very  large 

*  In  some  cases,  If  the  exponential  Is  approximated  by  rational  functions 
the  method  does  not  preserve  the  Lie  group  structure, 

in  particular,  when  Q  =  SL(n) 

==^  Another  class  of  methods  is  required. 
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2  Solvability  +  splitting 


The  procedure 

For  the  linear  system 


Y'=A{t)Y,  r(0)=/, 


we  denote  Yq  =  Y ,Aq=A  and  suppose  that 

^o(0  ~  ^0+  (0  +^0_  (t), 

where 

^0+  €  Vn  is  strictly  upper-triangular 
Ao_  G  A„  is  weakly  lower-triangular. 


NEW  NUMERICAL  INTEGRATORS  BASED  ON  SOLVABILITY  AND  SPLITTING  -  35 


C ◄ ►DOX 


2  Solvability  +  splitting  (II) 


The  idea  is  to  write  the  solution  as  a  product  of  upper  and  lower 


triangular  matrices. 
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2  Solvability  +  splitting  (II) 


The  idea  is  to  write  the  solution  as  a  product  of  upper  and  lower 
triangular  matrices. 

More  specifically,  we  propose  the  following  factorization: 


such  that 

=  io(0)  =  / 
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2  Solvability  +  splitting  (II) 


The  idea  is  to  write  the  solution  as  a  product  of  upper  and  lower 
triangular  matrices. 

More  specifically,  we  propose  the  following  factorization: 


Yo{t)  =  Loit)Uoit)Yiit) 


such  that 


L'q  =  Ao_  {t)Lo, 

Observe  then  that  Lo{t)  can  be  obtained 


Lo(0)=/ 

by  quadratures  and  Lo{t)  e  A„. 
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2  Solvability  +  splitting  (III) 


Now  we  form  the  matrix 


^0  '^0+-^ 
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2  Solvability  +  splitting  (III) 


Now  we  form  the  matrix 


^0  '^0+-^ 


which  can  also  be  split  as 

Co(r)  =  Co,(r)+Co_W, 


where 

Co+  €  Vn  is  weakly  upper-triangular 
Co-  €  A„  is  strictly  lower-triangular. 
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2  Solvability  +  splitting  (IV) 


Next  we  choose  Uq  as  the  solution  of 


Ul)  =  Co^{t)Uo,  Uo(0)=l 
so  that  t7o(t)  can  also  be  obtained  by  quadratures. 
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2  Solvability  +  splitting  (IV) 


Next  we  choose  Uq  as  the  solution  of 


t/i  =  Co+(f)t/o,  t/o(0)=/ 

so  that  C/o(0  can  also  be  obtained  by  quadratures. 

It  is  easy  to  show  that  Y\  satisfies 


Y[=Ax{t)Yx,  71  (0)=/, 


with 

A,  =C/„-‘Co_C/o- 
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2  Solvability  +  splitting  (V) 


This  gives  a  single  step  of  the  solvable  cycle,  which  we  repeat  with  Ai. 


Ai  —  Ai_|_ +Ai_,  Ai^  €  Vn5  €  A„ 


Yi=LiUiY2 

L[=AiM.  A(0)=/ 


etc. 
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2  Solvability  +  splitting  (VI) 


In  this  way  one  has  the  following  algorithm: 


Y  =  Yo=LoUQLxU^---LkUkYk+i 

with  (^  =  0, 1,2, . . .) 

■^k+  €  Vn?  -^k-  €  A 

L[=At_Lt,  Lt(0)=/ 

Q  =  +  Q_ 

Ct+  €  Vri) 

Ul  =  Ct^Ut,  Ut(0)  =  I 
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2  Solvability  +  splitting  (VII) 


and  finally 


=  t4  ^^k-Uk:  ^k+l  =M+\^k+\ 
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2  Solvability  +  splitting  (VII) 


and  finally 

=  Uj^^Ck_Uk,  =Ak+iYk-\-i 

Usually  the  factorization  is  truncated  by  taking  Yk+i  =  I. 
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2  Solvability  +  splitting  (VII) 


and  finally 


=  t4  ^^k-Uk:  ^k+l  =M+\^k+\ 


Usually  the  factorization  is  truncated  by  taking  =  I. 

In  what  follows  we  will  analyse  the  main  features  of  this  procedure  as  a 
numerical  integrator. 
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2.1  Order  of  the  method 


Suppose  that  A{t)  =  e(ao  +  H - )  for  some  parameter  e  >  0. 
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2.1  Order  of  the  method 


Suppose  that  A(f)  =  £{ao  +  ait  +  a2t^  H — )  for  some  parameter  8  >  0. 
Then 


Aj_  =  (eai  +t{ea2  +  8^a3)  +  0{t^)) 

Aj^  =  (ePi  +  r(8(32  +  e^ps)  +  0{t^)) 

for  7  =  1,2,.. .,  so  that 

Lj(t)  =  /H - - — (/8)”^'"'"^ai  H - - — /”^'"'"^8"j'(8a2  +  £^0C3)  H — 

Uj(t)  =  /+^_(,e)'”.+i|3i  +  _!_rj+V;(E|32  +  £2|33)  +  ... 

nij  +  l  rrij  +  2 
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2.1  Order  of  the  method  (II) 


Furthermore, 


rij+i  =  rij  +  mj  +  l 

mj+\  =  nj  +  2mj  +  2  j  =  l,2, ... 
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2.1  Order  of  the  method  (II) 


«7+l 


Mj+i 


Hj  +  my  +  1 

Yij  2my  “h  2  y  =  1 , 2 , . . . 


• 

J 

my 

1 

1 

2 

2 

4 

7 

3 

12 

20 

4 

33 

54 

5 

88 

143 
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2.1  Order  of  the  method  (III) 


In  consequence, 
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2.1  Order  of  the  method  (III) 


In  consequence, 


(1)  This  algorithm  could  be  useful  for  problems  of  the  form 


Y'  =  (Bo+zBi)Y 


if  the  system  Y'  —  BqY  can  be  solved  exactly. 
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2.1  Order  of  the  method  (III) 


In  consequence, 

(1)  This  algorithm  could  be  useful  for  problems  of  the  form 


Y'  =  {Bo  +  zBi)Y 


if  the  system  Y'  =  BqY  can  be  solved  exactly. 
(2)  The  order  of  approximation  is... 
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2.1  Order  of  the  method  (IV) 


7o  ~  LqUq 
7o  ~  LqUqLx 
Yq  w  LqUoL\U\ 

7o  ~  LqUoL\U\L2 
7o  «  L0U0L1U1L2U2 

Yo  «  LQU0L1U1L2U2L3 
Yo  «  L0U0L1U1L2U2L3U3 


is  order  1 

2 

4 

7 

12 

20 

33 
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2.1  Order  of  the  method  (IV) 


7o  ~  LqUq 
7o  ~  LqUqLx 
Yq  w  LqUoL\U\ 

7o  ~  LqUoL\U\L2 
7o  «  L0U0L1U1L2U2 

Yo  «  LQU0L1U1L2U2L3 
Yo  «  L0U0L1U1L2U2L3U3 


is  order  1 

2 

4 

7 

12 

20 

33 


With  only  4  solvable  cycles  we  get  order  33! 
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2.1  Order  of  the  method  (IV) 


7o  ~  LqUq 
7o  ~  LqUqLx 
Yq  w  LqUoL\U\ 

>0  ~  L0U0L1U1L2 

Yo  «  L0U0L1U1L2U2 
Yo  «  LQU0L1U1L2U2L3 
Yo  «  L0U0L1U1L2U2L3U3 


is  order  1 

2 

4 

7 

12 

20 

33 


With  only  4  solvable  cycles  we  get  order  33! 

...if  we  can  compute  and  C4  up  to  this  order... 
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2.2  Questions 


Several  problems  involved 
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2.2  Questions 


Several  problems  involved 


(1)  Does  the  approximate  solution  evoive  in  the  Lie  group  if  A  is  in  the  Lie 
algebra,  i.e.,  is  it  a  Lie  group  method? 
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2.2  Questions 


Several  problems  involved 

(1)  Does  the  approximate  solution  evolve  in  the  Lie  group  if  A  is  in  the  Lie 
algebra,  i.e.,  is  it  a  Lie  group  method? 

(2)  Solve  explicitly  the  systems  =Ak_Lk  and  =  Q+f4 
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2.2  Questions 


Several  problems  involved 


(1)  Does  the  approximate  solution  evoive  in  the  Lie  group  if  A  is  in  the  Lie 
algebra,  i.e.,  is  it  a  Lie  group  method? 

(2)  Soive  explicitly  the  systems  =  Ak_Lk  and  U^.  = 

(3)  Approximate  efficiently  the  (multiple)  integrals  involved. 
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3  Practical  issues 


(1)  Preservation  of  the  Lie-group  structure 


If  A(r)  G  5l{n),  the  algorithm  provides  by  construction  approximations  to 
Y{t)  in  SL(n). 
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3  Practical  issues 


(1)  Preservation  of  the  Lie-group  structure 


If  A(r)  G  5l{n),  the  algorithm  provides  by  construction  approximations  to 
Y{t)  in  SL(n). 

Proof.  Ak  =  Ak^  +Ak_ ,  with  Ak^  G  Vn.  M.  G  A„.  In  fact  Ak_  belongs  to 
a  solvable  subalgebra  of  s[(n).  Therefore  the  solution  of 


L[=Ak_Lk.  A(0)=/ 


Lkit)  G  SL(n)  (in  fact,  a  solvable  subgroup  of). 
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3  Practical  issues 


(1)  Preservation  of  the  Lie-group  structure 


If  A(r)  G  s[(n),  the  algorithm  provides  by  construction  approximations  to 
Y{t)  in  SL(n). 

Proof.  Ak  =  Ak^  +Ak_ ,  with  Ak^  G  Vn.  G  A„.  In  fact  Ak_  belongs  to 
a  solvable  subalgebra  of  5i{n).  Therefore  the  solution  of 


L[=Ak_Lk,  A(0)=/ 

Lk{t)  G  SL(n)  (in  fact,  a  solvable  subgroup  of). 

T^^{Ak+)  =  0,  and  the  trace  is  invariant  under  similarity,  so  that 

Tr(Q)  =  Tr(4-'At^Lt)  =  Tr(AtJ  =  0  ^  Q  €  sl(n) 
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3  Practical  issues  (II) 


Next,  Ck  =  Ck+  +Ck_,  with  Ck+  G  Vn.  Ck.  G  A„  and  Uk,  solution  of 


U'k  =  Ck^Uk,  Uk{())=I 


belongs  to  SL(n).  Finally 

At+i  =  U^'Ct_Ut  S5[(n) 


and  the  process  is  repeated. 
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3  Practical  issues  (II) 


Next,  Ck  =  Ck+  +Ck_,  with  Ck+  G  Vn.  Ck.  G  A„  and  Uk,  solution  of 


U'k  =  Ck^Uk,  Uk{())=I 


belongs  to  SL(n).  Finally 

At+i=Ut'CtUt  e5[(«) 

and  the  process  is  repeated.  □ 

Other  properties  (i.e.,  orthogonality)  are  preserved  only  up  to  the  order 
of  the  method. 
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3  Practical  issues  (III) 


(2a)  Explicit  solution  of  =Ak_Lk 
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3  Practical  issues 


(2a)  Explicit  solution  of  =  Ak_Lk 

Consider  ^  =  0  and  denote  Ao(?)  =  {aij),  ij  =  1 

j<i 


JO 


Then  the  solution  of  Lq  =  Ao_  {t)LQ,  Lo(0)  =  /  is 


i  —  1 , . . . ,  n 


.  ,  Tl,  Lq  (f  )  —  (Z^/ y ) , 


(3) 


0 


)^/:7  (^1 )  j 

k=j 


• _  • _  -I  •  -I 

I  ^  ^  ^  J  L  ^  ^  'I/  -Lb 
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3  Practical  issues  (IV) 


(2b)  Explicit  solution  of  Ul  = 


NEW  NUMERICAL  INTEGRATORS  BASED  ON  SOLVABILITY  AND  SPLITTING  -  69 


C ◄ ►DOX 


3  Practical  issues  (IV) 


(2b)  Explicit  solution  of 

Consider  ^  =  0  and  denote  Co(0  =  {cij),  ij  =  1 

j>i 

Cii{t)=  I  Cii{ti)dti. 
Jo 

Then  the  solution  of  Uq  =  Co^{t)Uo,  Uo{0)  =  /  is 
Unit)  = 


.  ,n,  Uo{t)  =  (Uij), 


(4) 


Uij{t)  =  e 


/  J 

Cu{t)  /  ^-Cuih) 


0 


(^1 )  (^1 )  j 

/:=/+! 


/  =  1 ,  .  .  .  ,  n  —  1 ,  J  =  /  +  1 ,  .  .  .  ,  n. 
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3  Practical  issues  (V) 


Explicit  expressions  for  the  elements  of  and  Uk  in  terms  of 


multivariate  integrals. 
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3  Practical  issues  (V) 


^  Explicit  expressions  for  the  elements  of  Lk  and  Uk  in  terms  of 
multivariate  integrals. 


They  can  be  evaluated  in  sequence  as  follows: 


1 


/  =  1 ,  .  .  .  , 

/  =  2 , . . . , 
/  =  3 , . . . , 


i  —  1 , . . . ,  /I 
i=  1 , . . . ,  n  —  1 

/  =  \  ^  ,  ^YL  —  2 


n 
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3  Practical  issues  (VI) 


In  principle,  the  integrals  appearing  in  and  C4  can  be  approximated  by 


quadrature  rules. 
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3  Practical  issues  (VI) 


In  principle,  the  integrals  appearing  in  and  C4  can  be  approximated  by 


quadrature  rules. 


To  minimise  the  computational  cost  this  has  to  be  done  by  using  the 
minimum  number  of  A  evaluations  in  each  integration  step. 
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3  Practical  issues  (VI) 


In  principle,  the  integrals  appearing  in  and  Uk  can  be  approximated  by 
quadrature  rules. 


To  minimise  the  computational  cost  this  has  to  be  done  by  using  the 
minimum  number  of  A  evaluations  in  each  integration  step. 


Question:  Is  it  possible  to  approximate  all  the  nested  integrals  with  the 
evaluations  required  to  compute 

Jo 


i.e.,  a  la  Magnus? 
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3  Practical  issues  (VI) 


In  principle,  the  integrals  appearing  in  and  Uk  can  be  approximated  by 
quadrature  rules. 


To  minimise  the  computational  cost  this  has  to  be  done  by  using  the 
minimum  number  of  A  evaluations  in  each  integration  step. 


Question:  Is  it  possible  to  approximate  all  the  nested  integrals  with  the 
evaluations  required  to  compute 

Jo 


i.e.,  a  la  Magnus? 


YES! 
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3.1  Example 


Illustration:  method  of  order  4  with  2  A  evaluations 


Step  /  =  0 1 — >  t  =  h. 
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3.1  Example 


Illustration:  method  of  order  4  with  2  A  evaluations 


Step  /  =  0 1 — >  t  =  h. 

1  -  Approximate  Aii{h),  i=  1, ...  up  to  order  4 


rh 

Aii{h)  =  /  aii{t)dt 

Jo 


—  2  (<3/;(0)  +Aaii{h/2)  +a/;(/r))  +  0{h^) 

=  Au{h)  +  0{h^) 


and  Aii{hl2),  i=  1, ...  ,n  —  1,  up  to  order  3  (necessary  to  approximate 

UjY 

Aii(h/2}  =  ^  (5a,, (0)  +  8a„(/>/2)  -  au(h}}  +  0(h*) 
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3.1  Example  (II) 


2-  Lii{h)  =  exp(A//(/i))  +  0(/z^)  (/  =  1, . . .  ,n)  and 
Lii{h/2)  =  exp{Aii{h/2))  +  0(/z"^)  (/  =  1, . . .  ,n  -  1). 


3-  Obtain  an  approximation  to  Lij{h),  j  <  i,  of  order  4  and  Lij{h/2)  of 
order  3 

Jo 


with 
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3.1  Example  (III) 


=  e-«('-)^(f;,(0)+4F,(V2)+F,W)  +  O(/.’) 


where  ^^7(0)  =  %(0)  and  Fij{h/2)  and  Fij{h)  have  to  be  approximated 
up  to  order  h^. 

The  sequence  of  computation  is  (/  =  2, . . .  ,n): 

(a)  Fi_i.,{hl2}  =  e-^«W2)aij_i{h/2)Li.ij.i{h/2)  +  0(h^) 

(b)  /v,,-i (A)  =  (A)  +  0{h^) 

(c)  Lij-i{h),  /  =  2, . . .  ,n  up  to  order  4 
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3.1  Example  (IV) 


L^-i(hl2)  =  ^»V>m 
(e)  Lij-2(h},  1  =  3,... 


^  (5a,-,_i  (0)  +  8F/,_i  (A/2)  -  F,/-i  (A))  +  0{h*) 
,n,  up  to  order  4  and  L,y_2(/r/2)  up  to  order  3 


...and  so  on. 
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3.1  Example  (IV) 


L,- i_  1  (h/l)  =  e^«<''/2) ^  ^ (o)  +  gir. (ft/2)  _  ,(*))  +  0{h^) 

(e)  Li^i-2{h),  i  =  3,...,  n,  up  to  order  4  and  L,y_2(/r/2)  up  to  order  3 


...and  so  on. 


In  this  way  we  have  Lo(/i)  computed  up  to  order  0{h^)  and  also  Lo{h/2) 
up  to  order  0{h'^)  with  2  evaluations  of  A(r). 
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3.1  Example  (V) 

3-  Next  we  compute  Cq: 


Co(0)  = 

Co{h/2)  = 

Co{h)  = 

Ao+(0)  error  0{h^) 

LQ\h/2)Ao^{h/2)Lo{h/2)  error  0(/z'^) 

{h)Ao^{h)Lo{h)  error  0{h^) 

4-  Cii{h)  —  I  (c;7(0)  +4cii{h/2)  +  0{h^) 

C,i(/!/2)  =  ^  (5cij(0)  +  Scu{h/2)  -  cuih))  +  0{h*) 
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3.1  Example  (VI) 


5-  i=  up  to  order  0{h^); 

Ui^i+i{h/2),  /  =  1, . . .  ,n  —  1,  up  to  order 
Uij+2{h),  i=  1, — 2,  up  to  order  0{h^); 
Ui^i+2{h),  i=  1, —  2,  up  to  order  0{h^)’, 


...  and  so  on. 

Thus  we  compute  ?7o(^)  with  error  0{h^)  and  also  Uo^h/l)  with  error 
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3.1  Example  (VII) 


6-Ai: 


Ai(0)  =  Co_(0)  error 

Ai(/z/2)  =  U^^{h/2)CQ_{h/2)UQ{h/2)  error  0(/z'^) 
Ai{h)  =  UQ\h)CQ_{h)Uo{h)  error  6>(/z^) 


...  and  the  process  is  repeated  again  for  the  second  cycle 
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3.1  Example  (VII) 


6-Ai: 


Ai(0)  =  Co_(0)  error 

Ai(/z/2)  =  U^^{h/2)CQ_{h/2)UQ{h/2)  error  0(/z'^) 
Ai{h)  =  UQ\h)CQ_{h)Uo{h)  error  6>(/z^) 


...  and  the  process  is  repeated  again  for  the  second  cycle 

it  is  possible  to  construct  a  method  of  order  4  with  only  2  A{t) 
evaluations  (3  for  the  first  step). 
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3.2  Other  possibilities 


One  could  use  other  quadrature  rules  instead,  for  instance 


Gauss-Legendre,  but... 
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3.2  Other  possibilities 


One  could  use  other  quadrature  rules  instead,  for  instance 


Gauss-Legendre,  but... 


in  that  case  there  are  not  enough  nodes  to  approximate  all  the 
multivariate  integrals  up  to  the  required  order. 
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3.2  Other  possibilities 


One  could  use  other  quadrature  rules  instead,  for  instance 
Gauss-Legendre,  but... 


in  that  case  there  are  not  enough  nodes  to  approximate  all  the 
multivariate  integrals  up  to  the  required  order. 

Remark.  This  is  not  the  case  with  Newton-Cotes  rules,  although  the 
error  introduced  may  be  important  for  higher  order  (negative  coefficients) 
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3.2  Other  possibilities 


One  could  use  other  quadrature  rules  instead,  for  instance 
Gauss-Legendre,  but... 


in  that  case  there  are  not  enough  nodes  to  approximate  all  the 
multivariate  integrals  up  to  the  required  order. 

Remark.  This  is  not  the  case  with  Newton-Cotes  rules,  although  the 
error  introduced  may  be  important  for  higher  order  (negative  coefficients) 

Solution:  use  G-L  with  matrix  evaluations  in  the  previous/next  step. 
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3.2  Other  possibilities 


One  could  use  other  quadrature  rules  instead,  for  instance 
Gauss-Legendre,  but... 


in  that  case  there  are  not  enough  nodes  to  approximate  all  the 
multivariate  integrals  up  to  the  required  order. 

Remark.  This  is  not  the  case  with  Newton-Cotes  rules,  although  the 
error  introduced  may  be  important  for  higher  order  (negative  coefficients) 

Solution:  use  G-L  with  matrix  evaluations  in  the  previous/next  step. 

method  of  order  4  with  2  evaluations  (and  1  in  the  next  step) 
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3.3  Some  methods 


Order  4 


y  «LoC/oLif/i 


*  Quadratures  NC  /  GL,  2  matrix  evaluations  per  step 

Order  6 

y  «Lot/oLit/iL2 


*  order  6  with  a  5  points  NC  quadrature  (4  evaluations  per  step) 

*  order  7  with  a  7  points  NC  (6  evaluations) 

Order  1 2 

Y  KLQU0UU1L2U2 


*  with  a  1 1  points  NC  (or  GL  involving  several  steps). 
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3.4  Variable  step  size 


Local  extrapolation  technique  is  trivial  to  implement  in  this  setting. 
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3.4  Variable  step  size 


Local  extrapolation  technique  is  trivial  to  implement  in  this  setting. 
For  instance, 


Fi  =  LoC/oLi 


Fi  =  LoUqLiUi=YiUi 


Then 

Yx-Yx=YxUi-Yi=Yi(Ui-I) 

and  II  —  Yi  ||  can  be  used  as  a  measure  of  the  error 
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3.5  Future  work 


*  Analyse  the  convergence  of  the  procedure 
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3.5  Future  work 
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*  Consider  numerical  examples  in  SL(n)  with  (very)  large  n 
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3.5  Future  work 


*  Analyse  the  convergence  of  the  procedure 


*  Consider  numerical  examples  in  SL(n)  with  (very)  large  n 


*  Highly  oscillatory  problems  (with  special  quadratures) 
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3.5  Future  work 


*  Analyse  the  convergence  of  the  procedure 

*  Consider  numerical  examples  in  SL(n)  with  (very)  large  n 

*  Highly  oscillatory  problems  (with  special  quadratures) 

*  Analyse  in  practice  the  preservation  of  other  structures  (Blanes  & 
Moan) 
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3.5  Future  work 


*  Analyse  the  convergence  of  the  procedure 

*  Consider  numerical  examples  in  SL(n)  with  (very)  large  n 

*  Highly  oscillatory  problems  (with  special  quadratures) 

*  Analyse  in  practice  the  preservation  of  other  structures  (Blanes  & 
Moan) 

*  Try  to  generalize  to  nonlinear  problems 
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The  End 
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